# Concurrent Sample Averages --------------------------------------------------------

historical_calibrations <- read_dta(paste0(EMISSIONS_DATA_IN, "/2. Historical Factors Preliminary Cleaned (Raw).dta"))

cems_lab_concurrent_scatter <-
historical_calibrations %>%
  # Only keep
  filter(phase=="II" & data_availability>.15) %>%
  mutate(pm_conc_cems_cal = case_when(
    device_type==2 ~ avg_cal_pm_cems,
    device_type==1 ~ (avg_cal_pm_cems * 1000000)/(normalised_flow * 3600)
  )) %>%
  filter(pm_conc_cems_cal>0 & pm_conc_cems_cal<2000 &
           pm_conc_lab_nm3>0 & pm_conc_lab_nm3<2000) %>%
  ggplot(aes(x=pm_conc_lab_nm3, y= pm_conc_cems_cal))+
  geom_point(color="dodgerblue", shape=1)+
  xlab("Sample PM (mg/Nm3)")+
  ylab("CEMS Reading (mg/Nm3)")+
  geom_smooth(aes(linetype="OLS"), method="lm", formula=y~x, se=F, color="black")+
  geom_abline(aes(linetype="y = x", slope=1, intercept=0),show.legend=F)+
  theme_classic()+
  scale_linetype_manual(values=c(1,3), name="")

ggsave(paste0(EMISSIONS_FIGS, "/Figure_F4.pdf"),
       cems_lab_concurrent_scatter,
       height = 4.5,
       width=6,
       units = "in"
       )
